
functions {
// scalar log-pdf of extreme observations other than last one
real GEVcpdf_lpdf(real x, real mu, real sig, real xsi){
		real z=(x-mu)/sig;
		real tau=exp(-6.0*abs(xsi));
		real lht;		
		if(1+xsi*z<tau){
			lht=-log(tau)/xsi-(z-(tau-1)/xsi)/tau;
		}
		else{
			lht=-log1p(xsi*z)/xsi;
		}
		return (1+xsi)*lht-log(sig);
	}
	
// scalar log-pdf of last (=smallest) extreme observation
real GEVpdf_lpdf(real x, real mu, real sig, real xsi){
		real z=(x-mu)/sig;
		real tau=exp(-6.0*abs(xsi));
		real lht;
		if(1+xsi*z<tau){
			lht=-log(tau)/xsi-(z-(tau-1)/xsi)/tau;
		}
		else{
			lht=-log1p(xsi*z)/xsi;
		}
		return (1+xsi)*lht-log(sig)-exp(lht);
	}

}	

data {
	real xi_min;
    real xi_max;
	int<lower=0> T;
	int<lower=0> nobs;
	matrix[nobs,T] y;      
}

parameters {

   // Level values for parameters
    real trans_xi_level;
    real ln_s_level;
    real m_level;
    
}

transformed parameters {
    
    
    //  Construct Basic Parameters
    real trans_xi = trans_xi_level;  
    real ln_alpha = 0.0;  // note that level of ln_alpha is fixed at 0
    real ln_s = ln_s_level;
    real m = m_level;

    // Transformed Parameters
    real alpha = exp(ln_alpha);
    real s = exp(ln_s);
   
    // Construct GEV Parameters
    real xi = xi_min + (xi_max-xi_min)*((exp(trans_xi))/(1+exp(trans_xi)));
    real sigma = s*(alpha^xi);    
    real mu = m + (s/xi)*((alpha^xi)-1);
}
	
model {

    trans_xi_level ~ normal(0.0,1.5);
    // Flat prior on ln_s_level and m_level

	for (t in 1:T) {
		for (j in 1:nobs-1) {
			y[j,t] ~ GEVcpdf(mu,sigma,xi);
		}
	y[nobs,t] ~ GEVpdf(mu,sigma,xi);
    }

}
